library(pacman)
p_load(
  here, fst, data.table, ggplot2, collapse, lubridate, yardstick, 
  patchwork
)
source('R/02-power-plants/predict-production.r')
theme_set(theme_minimal(base_size = 14))

# Function to propotionately split total excess load 
eload_proportinal = function(
  nerc_load_dt,
  quants = seq(0,1,by = 0.01),
  type = 'ic',
  manual_eload_limits = NULL
){
  if(type == 'ic'){
    # Total load max and min
    tot_eload_lim =
      nerc_load_dt[,
        .(tot_load = sum(excess_load)), 
        by = .(
          ic = fcase(
            nerc_adj %in% c('cal', 'wecc'), 'west',
            nerc_adj == 'tre', 'texas',
            default = 'east'
          ),
          datetime_utc
        )
      ][,.(
        max_eload = max(tot_load), 
        min_eload = min(tot_load)), 
        keyby = ic
      ]
    # Percent of total on average for each region
    nerc_adj_wts = 
      nerc_load_dt[,.(
        tot_region_load = sum(excess_load)), 
        keyby = .(
          ic = fcase(
            nerc_adj %in% c('cal', 'wecc'), 'west',
            nerc_adj == 'tre', 'texas',
            default = 'east'
          ),
          nerc_adj
        )
      ][,
        pct_load := tot_region_load/sum(tot_region_load), 
        by = ic
      ]
    # Scaling from min to max by the quants
    tot_eload_dt =
      tot_eload_lim[,
        as.list(data.table(quantile = quants)), 
        by = tot_eload_lim
      ][,tot_eload := min_eload + quantile*(max_eload-min_eload)]
    # Making load table 
    load_dt_raw = 
      merge(
        tot_eload_dt, 
        nerc_adj_wts,
        by = 'ic',
        allow.cartesian = TRUE
      )[,.(
        nerc_adj, 
        quantile,
        tot_eload,
        square = FALSE,
        excess_load = tot_eload*pct_load
      )]
  }else if(type == 'all'){
    # Total load max and min
    tot_eload_lim =
      nerc_load_dt[,
        .(tot_load = sum(excess_load)), 
        by = .(
          ic = fcase(
            nerc_adj %in% c('cal', 'wecc'), 'west',
            nerc_adj == 'tre', 'texas',
            default = 'east'
          ),
          datetime_utc
        )
      ][,.(
        max_eload = max(tot_load), 
        min_eload = min(tot_load)
      )]
    # Manually changing limits if given
    if(!is.null(manual_eload_limits)){
      tot_eload_lim[,':='(
        min_eload = manual_eload_limits[1],
        max_eload = manual_eload_limits[2]
      )]
    }
    # Percent of total on average for each region
    nerc_adj_wts = 
      nerc_load_dt[,.(
        tot_region_load = sum(excess_load)), 
        keyby = nerc_adj
      ][,pct_load := tot_region_load/sum(tot_region_load)]
    # Scaling from min to max by the quants
    tot_eload_dt =
      tot_eload_lim[,
        as.list(data.table(quantile = quants)), 
        by = tot_eload_lim
      ][,tot_eload := min_eload + quantile*(max_eload-min_eload)]
    # Making load table 
    load_dt_raw = 
      tot_eload_dt[, 
        as.list(nerc_adj_wts),
        by = tot_eload_dt,
      ][,.(
        nerc_adj, 
        quantile,
        tot_eload,
        square = FALSE,
        excess_load = tot_eload*pct_load
      )]
  }else if(type == 'quants'){
    load_dt_raw = 
      map_dfr(
        unique(nerc_load_dt$nerc_adj),
        \(nerc){
          data.table(
            nerc_adj = nerc, 
            quantile = quants,
            square = FALSE,
            excess_load = quantile(
              nerc_load_dt[nerc_adj == nerc]$excess_load, 
              probs = quants
            )
          )
        }
      )
    load_dt_raw[,tot_eload := sum(excess_load), by = quantile]
  }else{
    return("Incorrect type")
  }
  # Deviating each region by 1mwh at every quantile 
  load_dt_adj = 
    load_dt_raw[,
      as.list(data.table(deviation = c('baseline',unique(nerc_load_dt$nerc_adj)))),
      by = load_dt_raw
    ][,excess_load := fcase(
      deviation == nerc_adj, excess_load + 1,
      deviation == 'baseline' | deviation != nerc_adj, excess_load,
      default = NA
    )]
  # Adding squared term
  load_dt = 
    rbind(
      load_dt_adj,
      load_dt_adj[,.(
        nerc_adj, 
        quantile, 
        tot_eload, 
        square = TRUE, 
        excess_load = excess_load^2,
        deviation
      )]
    )
  return(load_dt)
}


# First making a table of loads we want to predict from -----------------------
  all_load_dt = 
    read.fst(
      path = here(
        "Data/electricity-generation/clean-load/nerc-load-dt-oge.fst"
      ),
      as.data.table = TRUE
    )[!(nerc_adj %in% c('AK','ASCC','HI')),
      .(nerc_adj = str_to_lower(nerc_adj), datetime_utc, excess_load)
    ]
  all_load_dt[,    
    ic := fcase(
      nerc_adj == 'tre', 'Texas',
      nerc_adj %in% c('wecc','cal'), 'West',
      default = 'East'
    ) |> factor(levels = c('West','Texas','East'))]
  quant_load_dt = 
    eload_proportinal(all_load_dt, type = 'quants')[,
      id := paste(deviation, ' - ', quantile), 
    ] 
  load_dt = prep_load_table(quant_load_dt[square == FALSE], id_var = 'id') #

# Step 2: Generate predictions from the model ---------------------------------
  set.seed(218)
  # Prepping data
  model_tables = prep_model_table()
  plant_dt = prep_plant_table(model_tables$model_dt)
  epsilon_dt = generate_epsilons(plant_dt)  
  md_path = 'Data/electricity-generation/plant-model-fit/marginal-damage'
  dir.create(here(md_path), showWarnings = FALSE)
  nerc_already_run = list.files(here(md_path)) |> 
    str_extract('cal|wecc|tre|serc|mro|npcc|rfc')|>
    str_to_upper() |>
    na.omit()
  # Running for each region 
  nerc_to_run = plant_dt[!(nerc_adj %in% nerc_already_run)]$nerc_adj |> unique()
  # Takes ~20 mins
  map(
    nerc_to_run,
    predict_production_region,
    load_dt = load_dt, 
    plant_dt = plant_dt, 
    model_tables = model_tables, 
    epsilon_dt = epsilon_dt, 
    path = md_path, 
    add_actual_data = FALSE, 
    id_var = 'id'
  )
  # Updating the already run list
  nerc_path = list.files(here(md_path), full.names = TRUE)
  nerc_already_run =  nerc_path |> 
    str_extract('(cal|wecc|tre|serc|mro|npcc|rfc)(?=\\.fst)')|>
    str_to_upper() |>
    na.omit()
  plant_gen_dt = 
    map(nerc_path, read.fst, as.data.table = TRUE) |>
    rbindlist()
  plant_gen_dt[,':='(
    deviation = str_extract(id, '.*(?=  -  )'),
    quantile = str_extract(id, '(?<=  -  ).*') |> as.numeric(), 
    id = NULL
  )]
  # Calculating marginal damages
  md_dt = 
    dcast(
      plant_gen_dt, 
      formula = plant_id_eia + quantile ~ deviation, 
      value.var = 'damages'
    )[,.(
      plant_id_eia, quantile, 
      md_cal = cal - baseline,
      md_mro = mro - baseline,
      md_npcc = npcc - baseline,
      md_rfc = rfc - baseline,
      md_serc = serc - baseline,
      md_tre = tre - baseline,
      md_wecc = wecc - baseline
    )]
  # Calculating marginal production
  mp_dt = 
    dcast(
      plant_gen_dt, 
      formula = plant_id_eia + quantile ~ deviation, 
      value.var = 'fit_net_generation_mwh'
    )[,.(
      plant_id_eia, quantile, 
      mp_cal = cal - baseline,
      mp_mro = mro - baseline,
      mp_npcc = npcc - baseline,
      mp_rfc = rfc - baseline,
      mp_serc = serc - baseline,
      mp_tre = tre - baseline,
      mp_wecc = wecc - baseline
    )]
  # adding back together 
  sim_plant_dt = 
    merge(
      plant_gen_dt[deviation == 'baseline',-'deviation'],
      md_dt, 
      by = c('plant_id_eia','quantile')
    ) |>
    merge(
      mp_dt, 
      by = c('plant_id_eia','quantile')
    )
  # Adding plant info and total load 
  sim_plant_dt =  
    merge(
      sim_plant_dt[,-c('fuel_category','nerc_adj','tot_eload')],
      plant_dt[,.(plant_id_eia, fuel_category, nerc_adj)], 
      by = 'plant_id_eia'
    ) |>
    merge(
      unique(quant_load_dt[,.(quantile, tot_eload)]), 
      by = 'quantile'
    )
  sim_plant_dt[,
    ic := 
      fcase(
        nerc_adj == 'TRE', 'Texas',
        nerc_adj %in% c("CAL",'WECC'), 'West',
        default = 'East'
      ) |> 
      factor(levels =c('West','Texas','East'))
  ] 
  # Summarizing by region and fuel type 
  region_sim_dt = 
    sim_plant_dt[,.(
      tot_gen_mwh = sum(fit_net_generation_mwh),
      tot_co2e_tons = sum(fit_emissions_tons_co2e),
      tot_nox_tons = sum(fit_emissions_tons_nox),
      tot_so2_tons = sum(fit_emissions_tons_so2),
      tot_pm25_tons = sum(fit_emissions_tons_pm25),
      tot_damages = sum(damages),
      co2e_tons_per_mwh = sum(fit_emissions_tons_co2e)/sum(fit_net_generation_mwh),
      nox_tons_per_mwh  = sum(fit_emissions_tons_nox)/sum(fit_net_generation_mwh),
      so2_tons_per_mwh  = sum(fit_emissions_tons_so2)/sum(fit_net_generation_mwh),
      pm25_tons_per_mwh = sum(fit_emissions_tons_pm25)/sum(fit_net_generation_mwh),
      damages_per_mwh = sum(damages)/sum(fit_net_generation_mwh)), 
      keyby = .(nerc_adj, ic, quantile, tot_eload)
    ]
  ic_sim_dt = 
    sim_plant_dt[,.(
      tot_gen_mwh = sum(fit_net_generation_mwh),
      tot_co2e_tons = sum(fit_emissions_tons_co2e),
      tot_nox_tons = sum(fit_emissions_tons_nox),
      tot_so2_tons = sum(fit_emissions_tons_so2),
      tot_pm25_tons = sum(fit_emissions_tons_pm25),
      tot_damages = sum(damages),
      co2e_tons_per_mwh = sum(fit_emissions_tons_co2e)/sum(fit_net_generation_mwh),
      nox_tons_per_mwh  = sum(fit_emissions_tons_nox)/sum(fit_net_generation_mwh),
      so2_tons_per_mwh  = sum(fit_emissions_tons_so2)/sum(fit_net_generation_mwh),
      pm25_tons_per_mwh = sum(fit_emissions_tons_pm25)/sum(fit_net_generation_mwh),
      damages_per_mwh = sum(damages)/sum(fit_net_generation_mwh)), 
      keyby = .(
        ic, 
        quantile,
        tot_eload
      )
    ]
  # Calculating marginal damages
  md_sim_dt = 
    sim_plant_dt[,
      lapply(.SD, sum),
      keyby = .(quantile, tot_eload),
      .SDcols = str_subset(colnames(sim_plant_dt),'md|mp')
    ] |> melt(
      id.vars = c('quantile','tot_eload')
    ) %>%
    .[,.(
      nerc_adj = toupper(str_remove(variable, 'm(d|p)_')),
      variable = str_extract(variable, 'm(d|p)'),
      ic = fcase(
          toupper(str_remove(variable, 'm(d|p)_')) == 'TRE', 'Texas',
          toupper(str_remove(variable, 'm(d|p)_')) %in% c("CAL",'WECC'), 'West',
          default = 'East'
        ) |> factor(levels =c('West','Texas','East')),
      quantile, 
      tot_eload,
      value
    )] |>
    dcast(
      formula = nerc_adj + ic + quantile + tot_eload ~ variable, 
      value.var = 'value'
    ) |>
    merge(
      quant_load_dt[
        square == FALSE & deviation == 'baseline', 
        .(tot_eload_ic = sum(excess_load)),
        by =.(
          quantile, 
          ic = fcase(
              toupper(nerc_adj) == 'TRE', 'Texas',
              toupper(nerc_adj) %in% c("CAL",'WECC'), 'West',
              default = 'East'
            ) |> factor(levels =c('West','Texas','East'))
          )
      ],
      by = c('ic','quantile')
    ) %>%
    .[,':='(
      nerc_adj_f = factor(
        nerc_adj,
        levels = c('CAL','WECC','TRE','MRO','NPCC','RFC','SERC')
      )
    )]

# Plots for the paper ---------------------------------------------------------
  region_md_demand_p = 
    ggplot(
      md_sim_dt[quantile %between% c(0.00,1.00)], 
      aes(
        x = tot_eload_ic/1e3, y = md/mp, 
        color = nerc_adj_f, linetype = nerc_adj_f
      )
    ) + 
    geom_line(linewidth = 1.25) +
    scale_color_brewer(name = 'Region',palette = 'Dark2') + 
    scale_linetype_manual(
      name = 'Region',
      values = c('solid','twodash','dashed','longdash','solid','dashed','solid')
    )+
    labs(
      x = 'Excess Load (GWh)',
      y = 'Marginal Damage ($/MWh)'
    ) +
    facet_grid(cols = vars(ic), scales = 'free_x')+ 
    theme(
      legend.position = 'bottom',
      legend.key.width = unit(3, 'cm')
    )
  region_md_demand_p
  ggsave(
    plot = region_md_demand_p,
    filename = here('figures/electricity-generation/eload-sim-region-md-demand.jpeg'),
    bg = 'white',
    width = 12, height = 5
  )
  region_mp_demand_p = 
    ggplot(
      md_sim_dt, 
      aes(
        x = tot_eload_ic/1e3, 
        y = mp, 
        color = nerc_adj_f, 
        linetype = nerc_adj_f
      )
    ) + 
    geom_line(linewidth = 1.25) +
    scale_color_brewer(name = 'Region',palette = 'Dark2') + 
    scale_linetype_manual(
      name = 'Region',
      values = c('solid','twodash','dashed','longdash','solid','dashed','solid')
    )+
    labs(
      x = 'Excess Load (GWh)',
      y = 'Marginal Production'
    ) +
    facet_grid(cols = vars(ic), scales = 'free_x')+ 
    theme(
      legend.position = 'bottom',
      legend.key.width = unit(3, 'cm')
    )
  region_mp_demand_p
  ggsave(
    plot = region_mp_demand_p,
    filename = here('figures/electricity-generation/eload-sim-region-mp-demand.jpeg'),
    bg = 'white',
    width = 12, height = 5
  )
  # Density of load
  eload_sim_ic_hist_p = 
    all_load_dt[,
      .(tot_eload_ic = sum(excess_load)/1e3),
      keyby = .(
        datetime_utc, 
        ic = fcase(
          nerc_adj == 'tre','Texas',
          nerc_adj %in% c('cal', 'wecc'), 'West',
          default = 'East'
        ) |> factor(levels = c('West','Texas','East'))
      )
    ] |>
    ggplot(aes(x =  tot_eload_ic)) + 
    geom_histogram(bins = 70) + 
    facet_wrap(~ic, scales = 'free_x') + 
    labs(
      x = 'Excess Load (GWh)',
      y = 'Number of Hours'
    )
  ggsave(
      plot = eload_sim_ic_hist_p,
      filename = here('figures/electricity-generation/eload-sim-ic-hist-p.jpeg'),
      bg = 'white',
      width = 12, height = 5
  )
  # Make 3 plots to get legend underneath 
  md_pW = 
    ggplot(
      md_sim_dt[quantile %between% c(0.0,1.00) & ic == 'West'], 
      aes(x = tot_eload_ic/1e3, y = md/mp, color = nerc_adj_f, linetype = nerc_adj_f)
    ) + 
    geom_line(size = 1.25) +
    scale_color_brewer(name = 'Region',palette = 'Dark2') + 
    scale_linetype_manual(
      name = 'Region',
      values = c('solid','twodash','dashed','longdash')
    ) +
    labs(
      x = '',
      y = 'Marginal Damage ($/MWh)',
      title = 'West'
    ) +
    ylim(75, 201) + 
    theme(
      legend.position = 'bottom',
      legend.key.width = unit(2.25, 'cm'),
      plot.title = element_text(
          hjust = 0.5,
          color = 'grey10',
          size = rel(1.0)
        )
    ) + 
    guides(
      color=guide_legend(nrow=2,byrow=TRUE),
      linetype=guide_legend(nrow=2,byrow=TRUE)
    )
    md_pT = 
      ggplot(
        md_sim_dt[quantile %between% c(0.0,1.00) & ic == 'Texas'], 
        aes(x = tot_eload_ic/1e3, y = md/mp, color = nerc_adj_f, linetype = nerc_adj_f)
      ) + 
      geom_line(size = 1.25) +
      scale_color_manual(
        name = '',
        values = RColorBrewer::brewer.pal(7, "Dark2")[3]
      ) +
      scale_linetype_manual(
        name = '',
        values = c('solid')#,'twodash','dashed','longdash')
      ) +
      labs(
        x = 'Excess Load (GWh)',
        title = 'Texas'
      ) +
      ylim(75, 201) + 
      theme(
        legend.position = 'bottom',
        legend.key.width = unit(2.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          color = 'grey10',
          size = rel(1.0)
        ),
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
      ) + 
    guides(
      color=guide_legend(nrow=2,byrow=TRUE),
      linetype=guide_legend(nrow=2,byrow=TRUE)
    )
    md_pE = 
      ggplot(
        md_sim_dt[quantile %between% c(0.0,1.00) & ic == 'East'], 
        aes(x = tot_eload_ic/1e3, y = md/mp, color = nerc_adj_f, linetype = nerc_adj_f)
      ) + 
      geom_line(size = 1.25) +
      scale_color_manual(
        name = '',
        values = RColorBrewer::brewer.pal(7, "Dark2")[4:7]
      ) + 
      scale_linetype_manual(
        name = '',
        values = c('solid','twodash','dashed','longdash')
      ) +
      labs(
        x = '',
        title = 'East'
      ) +
      ylim(75, 201) + 
      theme(
        legend.position = 'bottom',
        legend.key.width = unit(2.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          color = 'grey10',
          size = rel(1.0)
        ),
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
      )+ 
    guides(
      color=guide_legend(nrow=2,byrow=TRUE),
      linetype=guide_legend(nrow=2,byrow=TRUE)
    )
  # Adding all three plots together 
  ggsave(
      plot = md_pW + md_pT + md_pE,
      filename = here('figures/electricity-generation/eload-sim-region-md-mwh-demand.jpeg'),
      bg = 'white',
      width = 12, height = 5
  )
  
# How much is consumed within region vs exported 
md_region_sim_dt = 
  sim_plant_dt[,
    lapply(.SD, sum),
    keyby = .(quantile, tot_eload, nerc_adj),
    .SDcols = str_subset(colnames(sim_plant_dt),'md|mp')
  ] |> melt(
    id.vars = c('quantile','tot_eload','nerc_adj')
  ) %>%
  .[,.(
    nerc_adj_prod = nerc_adj,
    nerc_adj_demand = toupper(str_remove(variable, 'm(d|p)_')),
    variable = str_extract(variable, 'm(d|p)'),
    ic = fcase(
        toupper(str_remove(variable, 'm(d|p)_')) == 'TRE', 'Texas',
        toupper(str_remove(variable, 'm(d|p)_')) %in% c("CAL",'WECC'), 'West',
        default = 'East'
      ) |> factor(levels =c('West','Texas','East')),
    quantile, 
    tot_eload,
    value
  )] |>
  dcast(
    formula = nerc_adj_prod + nerc_adj_demand + ic + quantile + tot_eload ~ variable, 
    value.var = 'value'
  ) |>
  merge(
    quant_load_dt[
      square == FALSE & deviation == 'baseline', 
      .(tot_eload_ic = sum(excess_load)),
      by =.(
        quantile, 
        ic = fcase(
            toupper(nerc_adj) == 'TRE', 'Texas',
            toupper(nerc_adj) %in% c("CAL",'WECC'), 'West',
            default = 'East'
          ) |> factor(levels =c('West','Texas','East'))
        )
    ],
    by = c('ic','quantile')
  ) %>%
  .[,':='(
    nerc_adj_prod = factor(
      nerc_adj_prod,
      levels = c('CAL','WECC','TRE','MRO','NPCC','RFC','SERC')
    ),
    nerc_adj_demand = factor(
      nerc_adj_demand,
      levels = c('CAL','WECC','TRE','MRO','NPCC','RFC','SERC')
    )
  )] %>%
  .[,share_of_mp := mp/sum(mp), by=.(quantile,nerc_adj_demand)]

region_md_demand_own_p = 
  ggplot(
    md_region_sim_dt[
      nerc_adj_prod == nerc_adj_demand
     & quantile %between% c(0.00,1.00)
    ], 
    aes(
      x = tot_eload_ic/1e3, y = md/mp, 
      color = nerc_adj_demand, linetype = nerc_adj_demand
    )
  ) + 
  geom_line(size = 1.25) +
    scale_color_brewer(name = 'Region',palette = 'Dark2') + 
    scale_linetype_manual(
      name = 'Region',
      values = c('solid','twodash','dashed','longdash','solid','dashed','solid')
    )+
    labs(
      x = 'Excess Load (GWh)',
      y = 'Marginal Damage ($/MWh)'
    ) +
    facet_grid(cols = vars(ic), scales = 'free_x')+ 
    theme(
      legend.position = 'bottom',
      legend.key.width = unit(3, 'cm')
    )
  region_md_demand_own_p

# PLOT FOR PAPER: DEMAND FILLED BUT ONLY IN OWN REGION
  md_pW_own = 
    ggplot(
      md_region_sim_dt[nerc_adj_prod == nerc_adj_demand & ic == 'West'], 
      aes(x = tot_eload_ic/1e3, y = md/mp, color = nerc_adj_demand, linetype = nerc_adj_demand)
    ) + 
    geom_line(size = 1.25) +
    scale_color_brewer(name = 'Region',palette = 'Dark2') + 
    scale_linetype_manual(
      name = 'Region',
      values = c('solid','twodash','dashed','longdash')
    ) +
    labs(
      x = '',
      y = 'Marginal Damage ($/MWh)',
      title = 'West'
    ) +
    ylim(50,215) + 
    theme(
      legend.position = 'bottom',
      legend.key.width = unit(2.25, 'cm'),
      plot.title = element_text(
          hjust = 0.5,
          color = 'grey10',
          size = rel(1.0)
        )
    ) + 
    guides(
      color=guide_legend(nrow=2,byrow=TRUE),
      linetype=guide_legend(nrow=2,byrow=TRUE)
    )
    md_pT_own = 
      ggplot(
        md_region_sim_dt[nerc_adj_prod == nerc_adj_demand & ic == 'Texas'], 
        aes(x = tot_eload_ic/1e3, y = md/mp, color = nerc_adj_demand, linetype = nerc_adj_demand)
      ) + 
      geom_line(size = 1.25) +
      scale_color_manual(
        name = '',
        values = RColorBrewer::brewer.pal(7, "Dark2")[3]
      ) +
      scale_linetype_manual(
        name = '',
        values = c('solid','twodash','dashed','longdash')
      ) +
      labs(
        x = 'Excess Load (GWh)',
        title = 'Texas'
      ) +
      ylim(50,215) + 
      theme(
        legend.position = 'bottom',
        legend.key.width = unit(2.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          color = 'grey10',
          size = rel(1.0)
        ),
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
      ) + 
    guides(
      color=guide_legend(nrow=2,byrow=TRUE),
      linetype=guide_legend(nrow=2,byrow=TRUE)
    )
    md_pE_own = 
      ggplot(
        md_region_sim_dt[nerc_adj_prod == nerc_adj_demand & ic == 'East'], 
        aes(x = tot_eload_ic/1e3, y = md/mp, color = nerc_adj_demand, linetype = nerc_adj_demand)
      ) + 
      geom_line(size = 1.25) +
      scale_color_manual(
        name = '',
        values = RColorBrewer::brewer.pal(7, "Dark2")[4:7]
      ) + 
      scale_linetype_manual(
        name = '',
        values = c('solid','twodash','dashed','longdash')
      ) +
      labs(
        x = '',
        title = 'East'
      ) +
      ylim(50,215) + 
      theme(
        legend.position = 'bottom',
        legend.key.width = unit(2.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          color = 'grey10',
          size = rel(1.0)
        ),
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
      )+ 
    guides(
      color=guide_legend(nrow=2,byrow=TRUE),
      linetype=guide_legend(nrow=2,byrow=TRUE)
    )
  # Adding all three plots together 
  ggsave(
      plot = md_pW_own + md_pT_own + md_pE_own,
      filename = here('figures/electricity-generation/eload-sim-region-md-mwh-own-demand.jpeg'),
      bg = 'white',
      width = 12, height = 5
  )